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Abstract 


A fully two dimensional magnetohydrodyanamics code has been developed to predict self-field, steady-state 
MPD thruster performance. The governing equations are outlined, and methods of solution are presented. 
Model predictions are compared with experimental data for two thruster geometries over a range of discharge 
currents and mass flow rates. Model limitations are evaluated, and issues concerning quasi-steady versus 
steady-state thruster comparisons are discussed. 


Nomenclature 


an 

finite difference coeffs. 

V 

b 

thrust coefficient (N/A 2 ) 

w 

B 

magnetic field (T) 

a 

e 

electric charge (C) 


E 

electric field (V/m) 

In 

F 

thrust (N) 

r n 

A m 

j 

current density (A/m 2 ) 

1 

J 

discharge current (A) 

T]r 

ku 

Boltzmann const. (J/K) 

*ss 

m 

electron mass (kg) 

Kt 

M 

ion mass (kg) 

A 

m 

mass flow rate (kg/s) 

InA 

n 

number density (m~ 3 ) 

Mo 

P 

pressure (Pa) 

i 

P 

power (W) 

p 

r a 

anode radius (m) 

<7 

r c 

cathode radius (m) 

r e 

R 

gets constant (J/K-mol) 

<t> 

Rm 

magnetic Reynolds number 


3? 

mass flow ratio 


S 

source term 


T 

temperature (K) 

LJ 

V 

velocity (m/s) 

Ucc 

V ez 

exhaust velocity (m/s) 

n 


discharge voltage (V) 
weighting factor 
ionization fraction (= 1) 
finite difference coeffs. 
differential equation coeffs. 
differential equation coeffs. 
viscosity (kg/m-s) 
thruster efficiency 
frozen flow efficiency 
thermal conductivity (Watt/m-K) 
second coeff. of viscosity (kg/m-s) 
Coulomb logarithm 
permeability of free space 
density function 
mass density (kg/m 3 ) 
electrical conductivity (mho/m) 
electron collision time (s) 
electric potential (V) 
viscous dissipation function 
magnetic field streamline 
viscous force vector 
convergence factor 
electron cyclotron frequency (s^ 1 ) 
Hall parameter 


I. Introduction 


The magnetoplasmadynamic (MPD) thruster is an attractive candidate for orbit raising applications 1 as 
well as the high power, long duration missions envisioned by the national Space Exploration Initiative 2 . In 
its basic form, the MPD thruster consists of a cylindrical cathode surrounded by a concentric anode (Figure 
1). An arc struck between the electrodes ionizes a gaseous propellant, and the interaction of the current 
with self-induced magnetic fields accelerates the plasma to produce thrust. Steady-state MPD thrusters have 
been operated at power levels approaching 600 kW, while pulsed, quasi-steady devices have operated in the 
megawatt range 3 . The engine is robust, and designed to provide low, continuous thrust at specific impulse 
values between 1,000 and 10,000 s. 

MPD thruster performance is currently limited by low thrust efficiency in the operating regimes of interest 3 . 
Applied magnetic fields 3,4 and flared electrodes 5 have been shown to Influence MPD thruster efficiency 
and specific impulse, but for reasons poorly understood at present. In addition, severe electrode erosion is 
observed at high values of J 2 /rh, corresponding to the onset of voltage oscillations and unsteady thruster 
operation. Several mechanisms have been proposed to explain the onset of oscillations, including anode mass 
starvation 6 , flow choking due to enhanced back-EMF 7,8 , and the triggering of electrothermal and gradient 
instabilities as the plasma approaches full ionization 9-11 . Different operating conditions may trigger one or 
a combination of the proposed mechanisms, limiting thruster performance and lifetime. 

Although MPD thrusters have undergone extensive experimental development since the early 1960’s 312 , a 
comprehensive theoretical analysis has been slowed by the complex nature of the coupled electromagnetic and 
gasdynamic acceleration processes. Transient and steady-state numerical modeling using 1-D, quasi- 1-D, and 
quasi-2-D approximations of the fluid magnetohydrodynamic equations have provided valuable insights into 
self-field MPD thruster operation 13- 27 , but they are by definition constrained in their ability to predict global 
thruster performance. The removal of such limitations via fully 2-D and quasi-3-D modeling has recently 
become practical with the emergence of high speed computational facilities, allowing model validation and 
refinement using the existing experimental data base while in return establishing a theoretical basis to guide 
further experimentation;" 

A fully two dimensional, time dependent code for modeling self- field coaxial MPD thrusters is currently 
under development at the University of Stuttgart 28-30 . The model employs adaptive grid techniques for 
enhanced spatial resolution near electrode surfaces and incorporates important second order effects such 
as Hall currents, viscosity, and thermal conductivity. The plasma is assumed to be a fully ionized single 
fluid, but separate electron and ion temperatures are calculated to provide more accurate estimates of the 
thermal transport coefficients. The equations are cast into conservative finite volume form and iterated until 
a steady state condition is reached. Results compare favorably with experimental tests of the Stuttgart ZT-1 
and ZT-2 cylindrical thrusters at low values of J 2 /m, However, the code and has not been fully tested at 
values of J 2 jrn which correspond to predominantly electromagnetic regimes of thruster operation 31 . Because 
the code is being developed specifically to support the Stuttgart self-field plasma accelerator activities, no 
provision has been included for applied magnetic field effects 32 . 

Applied magnetic field effects are included in a cylindrical steady-state thruster code developed by Tanaka 
and Kimura 33 . The electromagnetic field equations are fully two dimensional, assuming axial symmetry for 
both the current distribution and the imposed applied magnetic field. The equations describing the plasma 
flow are based on a single-fluid, quasi- 1-D approximation, however, and viscous effects are neglected. In 
addition, the Hall parameter and electrical conductivity are assumed constant throughout the simulation 
region. The equations were written in a finite element formulation, and a Newton-Raphson iterative scheme 
was used to solve for the induced magnetic fields and current densities. Runge-Kutta algorithms were used to 
evaluate the one dimensional flow equations, and the entire process was repeated until the induced magnetic 
field and the plasma flow velocity coincided with their previous loop values. Calculations were performed at 
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discharge currents of 1000 A and 2000 A with a constant mass flow rate of 100 mg/s and applied magnetic 
field values up to 0.2 T. Although not compared with specific experimental results, the model predicted 
trends which had been observed experimentally 33 . 

A model which incorporates both self-field and applied magnetic field effects is currently being developed 
at the NASA Lewis Research Center to support a wide range of experimental MPD thruster activities. The 
fully 2-D steady-state code is based on the single fluid, single temperature MHD equations 3 , which assume 
axially symmetry but include azimuthal components of momentum, current density, and magnetic field. 
Viscous effects and Hall terms are included, as are non-constant expressions for the electrical conductivity 
and thermal conductivity. The plasma is currently assumed to be a fully (singly) ionized perfect gas.^ The 
governing MHD equations are written in finite difference formulation and iteratively solved on a uniform 
rectangular grid. The code is built using modular subroutines to allow for model expansion and refinement 
with minimal algorithm reconstruction. 

The self-field version of the NASA LeRC code is operational, and is discussed in Section II. Applied magnetic 
field terms, which do not contribute to self-field thruster operation, are retained in the derivative equations 
for completeness. Experimental comparisons with a coaxial, steady-state thruster and a flared, quasi-steady 
thruster are presented in Section III, and the paper concludes with a brief summary of results in Section IV. 


II. MPD Thruster Model 


The MPD thruster model presented below is based upon a single fluid, single temperature approxima- 
tion of the steady-state magnetohydrodynamic (MHD) equations, written in cylindrical coordinates with 
assumed symmetry about the centerline. Azimuthal components are retained, rendering the code fully two- 
dimensional. The plasma is currently assumed to be fully, singly ionized, and is described by a perfect gas 
equation of state. 

The fluid description is appropriate for the density and temperature regimes of interest in MPD thruster 
operation 35 , although such a formulation precludes a detailed evaluation of electrode fall potentials, and 
does not allow for potentially important kinetic effects such as plasma microinstabilities 30 . These limitations 
can be overcome by incorporating an electrode sheath model to provide boundary conditions for the fluid 
code (which in turn provides boundary conditions for the sheath model), and by enveloping the effects of 
microturbulence within modified expressions for the plasma transport coefficients. It is anticipated that such 
refinements will be included as appropriate models become available. 

The steady-state model is adequate for describing thruster behavior in regimes other than the brief transient 
start-up period or the onset of plasma oscillations associated with thruster instabilities. The startup period 
is, with proper design, a negligible fraction of the thruster lifetime and may be safely ignored in evaluating 
steady-state thruster performance. To properly model the onset of thruster instabilities as they evolve from 
steady-state operation would require a time-dependent plasma kinetics model and significant computational 
resources. An alternate and efficient near-term method is to develop separate models which predict regimes of 
operation under which such instabilities might occur, and incorporate those predictions as limiting conditions 
in the fluid code. MHD codes in turn can provide information on plasma conditions preceding the onset of 
thruster instability, as starting conditions for the instability models. While no single code can hope to model 
every effect, a proper combination of models can provide the understanding necessary for successful MPD 
thruster development. 

With this brief introduction, the equations incorporated into the self-field MPD thruster model are presented 
below. 

Ila. Electromagnetic Field Equations. 
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The basic set of electromagnetic equations includes the full complement of Maxwell’s equations: 


(a) V • B = 0 
(c) V E = p t /e 0 


0 


(b) V x B = poj 

(d) = . 


( 1 ) 


which incorporate the steady-state, single fluid plasma approximations. A general Ohm’s law of the form: 

n /r 


j = a E - x #)] — — (/ x s) 


( 2 ) 


is used to relate the current density j to the plasma velocity tf and the electric ( E ) and magnetic (B) fields. 
Ion slip terms are neglected due to the assumption of full ionization. The electrical conductivity a is given 
by the classical Spitzer-Harm conductivity for a fully ionized plasma 34 : 


cr = 1.53 x 10” 2 


rpZf2 

In A 


where T is the temperature in degrees-Keivin, InA is the Coulomb logarithm: 

lnA« 23 -ln' L22xW/2 ‘ 


T 3 / 2 


(3) 


(4) 


and n is the plasma number density expressed in particles per cubic meter. The Hall parameter Q is the 
product of the electron cyclotron frequency (w ce ) and the electron collision time (t c ): 


n = ui r . r. 


f eB \ / mer\ 


(5) 


= 9.6x10 


16 


/ t 3/2 b \ 

V^nlnA/ 


where m is the electron mass, e is the electron charge, and B is the magnitude of the local magnetic field. 
Equation 2 may be rearranged to solve for the electric field distribution: 



E = - 

1 

< 

.11 

are given 

by: 


E — 

d<t> 

1 

£j t — 

dr 

a 

Eq — 

84> 

~d8 

~ 0 

E x = 

d<t> 

dz 

1 

cr 


jr+ f0’ x *)]-(*■**) 


Jr + ^ (jeB, - j z B 0 ) j - (t > 9 B Z - v z 


n 

Jz + -jj ( Jr Bf) - j$B r ) 


Be) 


- (v r B e - v e B r ) (7c) 


( 6 ) 

(7a) 

(7b) 

(7c) 


The radial electric field is integrated from cathode to anode to find the potential drop <fi across the plasma. 
The divergence of Eq. 6 may be used to generate a Poisson equation for the potential, or assuming quasineu- 
trality a simple Laplace equation may be solved to give the potential distribution, with the calculated plasma 
fall providing the necessary boundary conditions. The first method is more rigorously correct, while the sec- 
ond method is computationally faster. Both methods were tried, and yield identical results for the plasma 
potential distribution. 
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A magnetic transport equation may be derived by combining Maxwell’s equations with the generalized Ohm s 
law: 


V x E = - 


dB 

dt 


„ Vx5 _ 

= V x + V x 

V / 


n 


Lmo crB 


((V x B) x £)j - V x (ifx b) =0 


Defining 0 — n/(<rB) and ip = tB 9 , the azimuthal component of the above equation can be written: 

d 2 ip dip dip . S ]2 ip c 

5^ +7l 57 +72 5 7 +73 ^ + ^ = 5 


( 8 ) 


( 9 ) 


where 


71 = 

72 = 

73 = - l Mo<7 


, 1 1 d(T crd/3 

- _ + t />-— +p 0 av r 

t a or r oz 


) 

dVr _ do«1"\ 

dr r dz \) 


(10) 


and the source term 5 is given by 


S — crrfi 


dB V 

dz 


a t 


(£♦$)( 
-f--S)(t 



( 11 ) 


For the self-field thruster, the source term reduces to zero. Once the magnetic transport equation is solved 
for ip, the value of B$ is readily obtained. The azimuthal magnetic field is a function of the total discharge 
current J, and through the relation ip = tBq ~ r(J/r) ~~ J, the function rp(r y z) =■ constant may be used to 
represent lines of total enclosed current. 


Returning to the Maxwell equations, the radial and axial current densities are obtained from the calculated 
magnetic field distributions: 


3 t 


it 


1 dB e 

Mo dz 



(12a) 

(12b) 


The azimuthal current density je vanishes due to the absence of applied radial or axial magnetic fields in 
the self-field MPD thruster. 


lib. Fluid Equations. 

The fluid equations are based on the conservation equations of mass, momentum, and energy. Conservation 
of mass for a compressible fluid is given by: 


V * {pv) = 0 


(13) 
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which can be written: 


djrpvr) djrpvz) _ ( 14 ) 

dr dz 

where p is the plasma mass density. Defining ( = rp, the mass conservation equation takes the simplified 
form: 

+ (15) 

dr dz 


which is solved for from which the value of p is determined. 

Momentum. The conservation of momentum, including viscosity, is expressed in vector form as. 

p(v-V)v = -Vp+ (Jx s) + ip (16) 


where p is the plasma pressure and $ is the viscous force vector: 

^=_-V[7 7 (V -ir)] + r 7 [V 2 ir+V(V tr)] + 2(V»j- V)v + Vrjx (V x v) 
3 

The viscosity rj for a fully ionized plasma is given by 3 : 


7] = 1.25 x 10 


- 19 _ 5 _ f MTukoT y' 2 f 2k B T \‘ 
8n'.4 \ 7r ) \ e 2 ) 


(17) 


(18) 


where M is the ion mass, m is the electron mass, e is the electron charge, k B is Boltzmann’s constant, n*_l 
for a singly ionized plasma, and, for a fully ionized gas (a = 1), 


A = 2 


In (l + or 2 ) - 


l + a 2 


= 0.386 


The radial component of the momentum equation is used to solve for the radial velocity v T : 


( dv T 
\ Vt dr 


+ v t 


dv T 

dz 



+ 


-TT + (j#*» - jM - v T 
or 

dv T dr] 4rj 3 2 v r d 2 v r 
J7dl + T dr 2 1 dz 2 


fin 1 dr, \ 4 9t>r fy drj\ 

1 3 r 2 + 3r dr ) 3 dr \ r dr ) 

7 ] d 2 v t dr) dv z 2 dr) dv t 
+ 3 drdz + dz~dr _ 3 dr dz 


(19) 


Equation (19) is recast in the form: 


d 2 v r 

77 


+r ^ +r ’ 


dv T 

dz 


- T 3 t v t + 


3rj d 2 v T 
4 dz 2 


where 


rl = 


r 


2 

r 


(V ,d]l _ 

\r + dr 4 ) 



-S T 


( 20 ) 


( 21 ) 
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and 


3 f dp . r] d 2 v x dv<hf_ _ i . ^9 \ 

Sr = 4 V Tr + liBl ~ 3xB * + 3 5rdz + dr 3 dr dz r ) 

The azimuthal component of the momentum equation is: 

f dv e , dv 9 ^v T v e \ _ _ . R 5^ / 5 *£ 

p(^v r — + ^— + — J - JxS r +77 +Jj 5z2 + ^ + 9r 

5t 7 dvs . Z' r] 1 

+ fca7 + ”*r^‘^J 

which yields an equation for u# of the. form: 

( 

*?- 

where 


■ 8Iv * + rJ^ + r|^-rS», + , 


dr 2 


dz 


d 2 ve 
dz 2 


= “^0 


r ; - (M--) 

r| - (g-~) 
rS - (*♦& + ?) 


and 


5*0 — jz Br ~ jr Bx 

Note that B T =B,=0 in the self- field thruster, and no azimuthal velocity will be generated. 

The axial component of momentum 

jj dr] \ dv^_ , 4 dr] dv, 

dz 


( 8v x dv,\ dp . . „ d 2 v, (n drA 5^ 4 5>7 

P { Vr ~d7 +V ‘~dI ) = -fc +JrBff ~ JffBr+V fr r+ (r + dr) dr + Z dz 


4 rj d 


T 3 dz 2 

yields an equation for the axial velocity v z : 


d 2 v , , .dv, , -,2 5V, 4q 5 2 u, 


dr 2 


+ r i^lli + r 2 

+ 1 * dr + z 


dz 3 dz 2 


= -S, 


where 


and 


« - (M--) 

^ ■ (II--) • 

S ' = -| +>B» - j.B, + ^ (^ + |) - I| ( V + dr ) 


, n d l v ,L 

+ 3 5r5z 


(22) 


(23) 


(24) 


(25) 


(26) 


V, 5Vr / 5tj\ _ 2 5»? / «r 5vr \ , V dV (27) 

z 2 + dz V3r + 5r / 3 dz\r dr ) 3 5rdz 


(28) 


(29) 


(30) 
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Energy. The conservation of energy, neglecting radiation losses, is given by: 

p (v *V) ( c p T ) ~ (v • V) P + V « (ktVT) + -f v x - j + <f> v 
where c p is the specific heat at constant pressure and kt is the Spitzer-Harm thermal conductivity 34 : 


kt ~ 4.4 x 10 


-10 


T 5/2 

In A 


The viscous dissipation is 34 ’ 35 : 


<t>v = V 2 


(£)’+(*) 


'+(ti 

1 dz 


+ 




dv e _ 
dr r 




where A represents the second coefficient of viscosity for a monatomic gas: 


A+j, = 0 


The energy equation is rewritten in the form: 

d 2 T L rl dT 


Krp 


dr 2 


where 


+ Tt dr 


r 1 — 

i r — 


F 2 — 

1 T — 


dT 

dz 


+ It — — h kt 0 , — —St 


d 2 T 

dz 2 


' kt dKT 

+ -5 pc p v, 

k r or 

'dK T \ 

jr-w) 


o 


and 


Sr = (v • V) P + -{- v x B^j * j 4* <f>v 

which is then solved for the single-fluid plasma temperature. 


(31) 


(32) 


^7 ( 33 ) 


( 34 ) 


(35) 


(36) 


(37) 


The set of fluid equations are closed with an ideal gas equation of state relating the plasma pressure, density, 
and temperature: 

p - pRT (38) 


l ie. Finite Difference Formulation 


The governing equations presented above have the general form: 


Or* 


(39) 


where the a n denote nonconstant coefficients and 5 represents a possible source term. The coupled form 
of the equations hinders an analytic solution, and suggests the use of an iterative numerical approach. The 
equations are written in finite difference form, and solved on a uniform rectangular grid using a generalized 
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Newton- Raphson iteration method 38 . Second derivatives are represented by a second-order central difference 
formulation: 

Y[j,i+l)-2Y(j, i) + Y{j,i-l) (40a) 

(Ar ) 2 

Y(j + l,i)-2Y(j,i) + Y{j - 1,*) ( 40b ) 

(A zf 


d 2 Y 

dr 2 

d 2 Y 

dz 2 


where the index j represents nodal points along the axial grid direction, i represents nodal points along the 
radial grid direction, and A r and A z are the radial and axial grid spacings, respectively. First derivatives 
are written using a first-order switching scheme, which is based upon the sign of the coefficients multiplying 
the derivatives 38 : 


dY 

dr 


dY_ 

dz 


At 

a 2 > 0 

(41a) 

At 

a 2 < 0 

(41b) 

YJj± 

Az 

Ot 3 > 0 

(41c) 

YU.i)-YU-U) 

<*3 < 0 

(41d) 

Az 



This scheme preserves the dominance of the Y(j,i) term and assures stability of the numerical solution. 
The first-order continuity equation (15) is solved for ( using an automatic upwind/downwind differencing 
method: 


dtv T 

dr 

d£v z 


dz 

where 


-i- [(tty - |uy|)£(j, i + 1) + {u f + |u/lK(i, i) + (-«» + M)£(j> 0 - («* + MW ’ * ~ < 42a - 1 

2Ar 


[(uy — |Uy|)£(j + l f t) + {vj + |v/|)^(j, i) + (— ' Vb + |Vft|)£(ji *) ( v b + |vfc|)^(i 1,0] ( 42b ) 


2Az 


Uf = 

^ (t>r(j,» + 1 ) + «r(ji0) 

(43a) 

Ufc = 

= (»rO’i i) + v T (j,i-l)) 

(43b) 

Vf = 

+ 1,0 + v *{j,i)) 

(43c) 

Vb = 

^ (Vr(j,i) + V T {j - 1 , 0 ) 

(43d) 


This method allows the continuity equation to be kept in conservative form while retaining numerical stability. 
Once written in their finite difference analogs, each equation is regrouped into the general form: 

F = do Y(j, i) + diF(j, 1+1)4- o 2 7 (j, i - 1) + a 3 Y(j + 1,0 + a^Y(j — 1, i) + S{j, i) (44) 

where the are nonconstant coefficients, which must be evaluated at each grid location. The generalized 
Newton- Raphson iteration method 

F 


Y (j, i) = Y{j,i)-u> 


[dF/dY{j,i)} 


(45) 
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is then used to iteratively solve each finite difference equation. An over-relaxation factor u is used to speed 
convergence; its value lies between 1 and 2, and was determined through trial and error for each equation 
being solved. 

Program execution is presented schematically in Figure 2. The fixed grid is divided into 50 radial by 100 
axial nodes, and each equation is iterated first in the axial and then in the radial direction until a relative 
convergence criteria is satisfied (typically to within a tolerance of 10” 4 ). The process is repeated for each 
equation in turn, and the entire loop through the full set of equations is repeated until the exhaust velocity 
and plasma potential have each converged to within 1% of their previous loop values. 

Mass flow conservation is checked by numerically integrating pv z radially along the thruster inlet and again 
at the thruster exit plane. Mass flow rates at the exit plane were initially found to be off by 15-20%, a 
consequence of the finite difference formulation employed; similar results have been reported for the Stuttgart 
finite difference MPD codes 28,30 . In an attempt to improve mass flow conservation, a weighting scheme was 
incorporated into the finite difference solution of the continuity equation. Equations (43) were replaced by: 


c 

Si 

1! 

ci 

-4 

h i + 1) + v r 

(j, i)™) /(I + tl>) 

(46a) 

Ub = (v T (j,i) + v r (j,i 

— l)u;) /(I 4- w) 

(46b) 

v ] = ( V *U + M) + v * 

(j, i)w)/{l + u>) 

(46c) 

Vb = {v T {. 

h i) + v r {j - 

l,i)w)/(l + ui ) 

(46d) 


where w is a weighting factor, initially set to unity. After each loop through the continuity equation, the 
ratio 3? of the mass flow at the thruster inlet plane to the mass flow at the exit plane was calculated and the 
weighting factor modified according to: 

w — 3?u> (47) 


The continuity equation was then reiterated until the mass flow ratio 3? = 1 ± 0.01. Although the value of 
the weighting factor changes over each loop, the actual velocities are not altered and the calculated mass 
flow rates can be brought into close agreement. 

Thrust and Specific Impulse. The electromagnetic thrust F em is calculated by numerically integrating the 
axial j x B body force over the current carrying volume, including regions downstream of the thruster exit: 

F em = [ (j x B) z dV (48) 

Jv 

Additional pressure forces are numerically integrated over the inner surfaces of the thrust chamber, and com- 
bined with the electromagnetic thrust calculation to provide total thrust. A second total thrust calculation 
is then performed using: 

F — 77lV ex -f* f ( j X B') z dV ou t ~h (Pe (^*0 

JVo»t 


The exhaust velocity u cx is given by: 


EiP0'a> i K(i a > 0 

Ex P( 0 


(50) 


where the axial index ja denotes the anode exit plane and the radial integration extends from the thruster 
centerline to the inner anode radius. The integration of the electromagnetic body force is performed over 
the current carrying volume external to the thruster. The pressure force term corresponds to an imbalance 
between the pressure at the anode exit plane and the background gas pressure, evaluated over the thruster 
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exit area A e . The complementary thrust calculation methods provide an independent check on the total 
thrust, and allow estimates to be made of the various thrust components. 

lid. Starting Values and Boundary Conditions 

Primary input into the code consists of the discharge current, propellant ion mass (amu), and propellant mass 
flow rate. Secondary input defines code options, such as maximum bulk electrode temperatures, background 
gas pressures, and initial estimates for the plasma transport coefficients. Bulk electrode temperatures are 
typically set at 3000 K, and electrode surface temperatures are then evaluated by averaging the temper- 
ature immediately above the electrode surface with the bulk electrode temperature, providing a smoother 
temperature transition between the electrode surface and the hot plasma. Background gas pressures are set 
between 1.3xl0~ 2 and 1.3xl0“ 4 Pa (10“ 4 and 10" 6 ton, respectively). No appreciable difference is found 
in calculated thruster performance using either background pressure value. 

The inlet velocity uq is assumed to be sonic and uniform at the backplate. The inlet density po is assumed 
to be uniform and is calculated from the 1-D continuity equation, m = PuUq^Oj where Ay is the exposed 
backplate surface area. Radial and azimuthal velocities are initially set to zero, and are defined by symmetry 
to remain zero along the centerline. Initial temperatures throughout the calculation region are set to the 
bulk electrode temperatures. Initial velocities and densities in the calculation region are set to their values 
at the backplate. A no-slip velocity condition is employed at all insulator and electrode surfaces. 

Electric fields at insulator surfaces are set equal to zero. Electric fields which enter perpendicular to electrode 
surfaces satisfy the condition dE n /dh—0, where h denotes a unit vector normal to the surface. Electric fields 
lying along electrode surfaces are set to zero, consistent with the electrodes being equipotential surfaces. The 
magnetic stream line rf) is set equal to -/i 0 J/(27r) along the backplate, zero along the centerline, and is assumed 
to be continuous at both the outer radial grid boundary and and downstream axial grid boundary. Setting 
\j )— 0 at the thruster exit plane prevents current from blowing out of the thruster, while setting 0=0 at the 
downstream grid boundary artificially compresses the current blown downstream. 

lie. Stability and Performance 

Steep radial pressure gradients, typically caused by a combination of high densities and sharp temperature 
gradients near the cathode tip, may cause the code to become unstable. Negative radial pressure gradients 
near the plume centerline can drive the radial velocity positive, causing the centerline density to plummet 
and the local Hall parameter to skyrocket. This effect can propagate into the field equations, precipitating 
numerical instabilities from which the steady-state code cannot recover. This problem is currently handled 
by under-relaxing the pressure gradient terms during the first several iterations, allowing the flow and 
discharge pattern to stabilize before the pressure gradient terms are brought up to their full values. A more 
elegant solution, currently under investigation, is to use nonuniform or multigrid algorithms to refine the 
computational mesh near the electrode surfaces, allowing the sharp gradients to be handled more effectively. 

To shorten loop convergence times, the temperature and velocity components are averaged with their previous 
loop values. This damps out minor fluctuations as a solution is approached, allowing more rapid loop 
convergence. Tests performed with and without smoothing yield identical results. Model convergence with 
smoothing is typically obtained within 10-12 loops through the equations, requiring approximately 45 CPU 
minutes on a VAX-9000 computer. The code has not yet been optimized for a CRAY computing system. 

III. Comparison with Experiment 


The self-field model was compared with published experimental results for the University of Stuttgart coaxial 
ZT-1 MPD thruster 39 and for the Princeton University half-scale flared anode thruster (HSFAT) 40 . The ZT- 


11 



1 thruster anode is segmented, allowing anode current distributions to be evaluated during steady state 
operation. In addition, the thruster has been numerically modeled by Sleziona ei a/. 28 , allowing direct 
comparisons between the Stuttgart and NASA LeRC codes. The HSFAT was chosen to test the accuracy 
with which a fixed numerical grid can model a flared thruster geometry, while the extensive experimental 
data which exists for the HSFAT 40 ’ 41 provides an excellent opportunity to exercise the code over a wide 
range of currents and mass flow rates. 

Ilia. ZT-1 Thruster Comparisons 

The University of Stuttgart ZT-1 MPD thruster is a cylindrical, steady-state, water-cooled research device 
with a 2% thoriated tungsten cathode and a segmented copper anode (Figure 3a). The two segments closest 
to the backplate are electrically isolated, with the remaining three rings connected to form the anode. The 
backplate is made of copper, and is electrically isolated from the electrodes. Propellant gas may be introduced 
through a backplate annulus surrounding either the cathode or the anode, or the flow may be split and sent 
through both injector rings. The anode inner diameter is 7.0 cm and has a total length of 15 cm, including 
the insulated anode segments. The cathode is 1.8 cm in diameter, and for the comparison below has an 
exposed length of 15 cm within the thrust chamber. 

To provide comparison data for the Stuttgart code, the ZT-1 thruster was operated at a discharge current 
of 6000 A, with an argon propellant mass flow rate of 6 g/s split 2:4 between anode and cathode 28 . Current 
distribution was measured in the three conducting anode segments, along with the total discharge voltage. 
The engine has not yet been mounted on a thrust stand, so no direct thrust measurements were available for 
comparison. Results of the Stuttgart model show excellent agreement with the current distribution along the 
anode, while the calculated discharge voltage, neglecting electrode fall potentials, is only slightly under the 
measured total voltage. The experimental and calculated values are summarized in Table 1, with predicted 
current contours displayed in Figure 3b. The thruster was operated steady-state, hence no direct current 
contour measurements could be performed. A similar thruster, designated ZP-1, was constructed to operate 
in a quasi-steady mode 39 ; comparisons between the measured current distributions and calculated current 
distributions were quite poor, however, and were ascribed to variations in the cathode emission distributions 
between quasi-steady and steady-state operation 30,39 . Unfortunately, the reported comparisons were at 
different discharge currents and mass flow rates (although similar / 2 /m), causing some ambiguity in their 
conclusions. 



U. STUTTGART 

NASA LERC 

ZT-1 THRUSTER 

EXPERIMENT 

MODEL 

MODEL 

ANODE SEGMENT # 

ANODE CURRENT FI 

R.ACTION 

1 

0.46 

0.44 

0.47 

2 

0.27 

0.27 

0.23 

3 

0.27 

0.29 

0.30 


DISCHARGE VOLTAGE (V) 

VOLTAGE (V) 

19 V 

14 V 

3.7 V 


Table 1: Current distribution and discharge voltage comparisons for the ZT-1 MPD thruster operating at 
6000 A, 6 g/s. Calculated voltages do not include electrode fall potentials. 


Results of the NASA LeRC simulation of the ZT-1 thruster are also shown in Table 1, and indicate fair 
agreement with both the experimental and numerical Stuttgart results. The steady-state LeRC code pre- 
dicts slightly higher current concentrations along the upstream and downstream anode sections, and a lower 
current concentration along the central anode section. The calculated enclosed current contours are dis- 
played in Figure 3c, and compare favorably with the calculated Stuttgart results. The different downstream 
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attachments shown in Figures 3b and 3c are a result of the boundary conditions used on the anode face. 
The Stuttgart model assumes a partially insulated anode face, which forces the current to reattach near the 
anode lip. In Figure 3c, the anode face is left uninsulated, and the downstream current does not appear to 
reattach within the simulated region. 

The plasma potential calculated with the LeRC code is substantially lower than the total measured voltage, 
reaching only 3.7 volts, or roughly 20% of the discharge voltage. The plasma drop, calculated by radially 
integrating Equation 7a from cathode to anode, does not include the electrode fall voltages or any plasma 
ionization penalties; however, neither does the Stuttgart calculation, which comes remarkably close to match- 
ing the full discharge voltage. The reason for this discrepancy lies in the different expressions used for the 
electrical conductivity. The average electrical conductivity calculated by the Stuttgart code is 29 2000-2500 
mho/m within the thruster, while the average value calculated in the NASA LeRC model is approximately 
6000 mho/m within the thruster. The calculated plasma potential is inversely proportional to the conduc- 
tivity (Equation 7a), which results in the higher potential calculated by the Stuttgart model. Resolving this 
issue will require better approximations for the plasma transport coefficients, as well as a sheath model to 
accurately calculate electrode fall contributions to the total voltage drop. 

Calculated thrust values were not reported for the ZT-1 thruster under the assumed operating conditions. 
A thrust of approximately 28 N, including 14 N of electromagnetic thrust, was calculated by the Stuttgart 
code 28 for the ZT-1 thruster using a discharge current of 10 kA and an argon mass flow rate of 6 g/s. The 
LeRC code was operated with these same parameters to verify the thrust calculation routine. The total 
calculated thrust was 27 N, with 13 N of electromagnetic thrust, in good agreement with the Stuttgart 
predictions. For the ZT-1 thruster operating at a discharge current of 6 kA, the LeRC code calculated a 
total thrust of approximately 11 N, with 4.7 N of electromagnetic thrust. The corresponding specific impulse 
is approximately 185 s. Using the measured discharge voltage of 19 volts at 6000 A and the calculated thrust 
of 11 N, the estimated efficiency of the ZT-1 thruster under the given operating conditions is approximately 
9%. 

Temperature, pressure, and density contours calculated with the NASA LeRC code are displayed for the ZT-1 
thruster in Figures 4(a) through 4(c). The pressure and density contours are normalized to their average inlet 
values and plotted as logi 0 [P/Po] and \ogi 0 [p/po], respectively, to bring out details in the exhaust plume. 
Electric equipotentials are plotted in Figure 4d, showing equipotential line convergence at the interface 
between the insulating and conducting anode segments, and at the inlet where the cathode penetrates the 
insulating backplate. Local Mach number contours are displayed in Figure 4e, and the velocity flow field 
is shown in Figure 4f. Of particular interest is the high temperature interelectrode region shown in Figure 
4a; a similar temperature increase was observed in the Stuttgart modeling results 28 ’ 29 . The convergence of 
equipotential lines along the junction of the insulated and conducting anode segments indicates a region of 
strong electric fields, and a plausible explanation for the increased temperature is localized Joule heating (cf 
Equation 31). To test this hypothesis, a second numerical run was performed for the ZT-1 thruster using 
identical operating conditions, but with the insulated segments replaced by conducting anode segments. 
Results are shown in Figures 5a through 5f, with the enclosed current contours displayed in Figure 5g. 
The high temperature interelectrode region appears in nearly the same location with approximately the 
same temperature as in the ZT-1 thruster with insulating segments. The high temperature region is thus 
not a result of localized Joule heating near the insulator/electrode interface, but is generic to the ZT-1 
thruster geometry under the assumed operating conditions. Further modeling would be useful to determine 
how different mass flow rates and discharge currents affect the interelectrode temperature distribution, and 
whether the localized increase in temperature can be reduced by changing the anode radius or the thruster 
length. A combination of experiment and simulation could determine the effect of the high temperature 
anomaly on electrode erosion, and produce design modifications which may mitigate the severe upstream 
cathode erosion observed in the ZT-1 thruster 42 . 

Temperatures are also observed to increase in the flow exhaust region, with maximum temperatures occuring 
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downstream of the thruster exit plane. The corresponding pressure distributions (Figures 4b and 5b) show 
increased pressures at locations of increased temperature, both within the chamber and downstream of the 
cathode tip. The pressure distribution for the ZT-1 thruster with insulating segments is somewhat more 
complex than the ZT-1 without insulators, indicating a pressure rise slightly upstream of the anode exit 
plane compared with a relatively smooth pressure transition in the ZT-1 without insulators. Mach number 
contours (Figures 4e and 5e) show the flow is more smoothly accelerated in the ZT-1 without insulators, 
although some flow deceleration occurs in the high pressure interelectrode region for both cases. Smoother 
flow is also shown in the velocity vector diagrams (Figures 4f and 5f), which display an interesting difference 
in the exhaust velocity profiles for the two cases. Calculated thrust values were 11.1 N for each case, including 
4.7 N of electromagnetic thrust. 

The ZT-1 thruster comparisons allow the self-field model to be tested against experimental and independent 
numerical results, and provide some confidence that the code can accurately predict major aspects of steady- 
state MPD thruster performance. However, the disparity in predicted plasma voltages underscores the 
need for accurate electrode sheath models and improved expressions for the plasma transport coefficients in 
calculating total discharge voltages. The prediction of a high temperature region occuring within the ZT-1 
thruster channel, and the possible correlation of this region with a measured increase in cathode erosion, 
reiterate the need for a balanced experimental and numerical approach in evaluating MPD thruster dynamics. 


Illb. HSFAT Thruster Comparisons 


The half-scale flared anode thruster is depicted schematically in Figure 6. Vacuum tank constraints limit 
thruster operation to pulsed, quasi-steady operation, and the continuous copper anode is not actively cooled. 
The 0.5 cm diameter cathode is made of 2% thoriated tungsten, and extends the full 10 cm length of the 
anode. The anode inner radius is approximately 1.6 cm at the backplate, extending 1.4 cm downstream 
before beginning a 12.5° half-angle flare. The flared section continues for roughly 4 cm, followed by a 4.6 
cm straight channel section. Propellant is injected through an insulating boron nitride backplate, with half 
the propellant mass flow injected through an annulus at the cathode base, and the remaining half injected 
through a ring of 12 evenly spaced ports radially located at 1.5 cm. Argon propellant was used in all tests, 
with mass flow rates of 0.75, 1.5, and 3 g/s and discharge currents up to 23 kA at the higher mass flow rate. 
Voltage-current characteristics were adequately described by the relation: 


V = 2.05 x 10' 3 J + 1.5 x 10“ 14 



(51) 


where V is the total discharge voltage in volts, J is the discharge current in amps, and m is the propellant 
mass flow rate expressed in kg/s. The onset of operational instability, defined as a 10% hash in the otherwise 
uniform voltage signal, was noted to occur at a critical J 2 jrn value of approximately 1.8 xlO 11 A 2 -s/kg. 


Table 2 presents a matrix of the modeled J 2 jrn values as a function of mass flow rate and total discharge 
current. Model comparisons were performed for each of the given mass flow rates at J 2 jrn values of approxi- 
mately 2.5 x 10 10 , 5.0 x 10 10 , L05 x 10 u , and 1.5 x 10 u A 2 -s/kg. Additional calculations were performed at 3 
g/s, J 2 /m = 6.1 x 10 9 A 2 -s/kg, and 1.5 g/s, J 2 jrn = 2.1 x 10 u A 2 -s/kg, to extend the range of comparisons. 

Thrust. Figure 7 presents measured and predicted thrust for the HSFAT as a function of current. The code 
accurately models the observed thrust for each of the mass flow rates and discharge currents up to J jrn 
values of approximately 1.5X10 11 A 2 -s/kg, at which point the calculated thrust begins to deviate from the 
observed values. As noted previously, the experimentally determined value at which 10% voltage oscillations 
occur is approximately 1.8 x 10 11 A 2 -s/kg, and it is probable that the simple single-fluid MHD code cannot 
accurately model the physical processes occuring in the thruster as operational instability is approached. It 
is interesting to note that the observed thrust at the 3 g/s mass flow rate displays a cubic dependence on the 
current above 1.8xl0 4 A, a trend which is missed by the code. These effects are a valuable aid in defining 
regimes of model accuracy, and provide essential benchmarks against which numerical simulations must be 
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HSFAT 


DISCHARGE CURRENT (xlO 3 A) 

MASS FLOW 

4.3 

6.1 

8.7 

10.0 

12.4 

15.0 

17.8 

21.3 

RATE (kg/s) 



J 2 jrn x 

L0 9 A 2 - 

s/kg 



3x 10 -3 

6.1 

- 

25 

- 

51 

- 

105 

150 

1.5x 10 -3 

- 

25 

50 

- 

103 

150 

212 

- 

7.5 x 10~ 4 

25 

50 

103 

150 

- 

- 

- 

- 


Table 2: Values of mass flow rate (m), discharge current {J), and J 2 /m (x 10 9 A 2 -s/kg) used in the HSFAT 
comparisons. 


tested and refined. 

Figure 8 displays the calculated electromagnetic thrust and the calculated total thrust plotted against J 2 jm 
for each of the propellant mass flow rates. The electromagnetic thrust contribution increases with increasing 
J 2 /rh t and dominates the acceleration at J 2 /m values above 8-10 xlO 10 A 2 -s/kg. At the highest J 2 /m 
value, the electromagnetic thrust calculated with Equation 48 is actually 2-5% higher than the total thrust 
calculated using Equation 49. Gilland 41 reports that higher thrust values were obtained for the HSFAT 
using electromagnetic probe results versus actual thrust measurements, and accounts for this difference by 
incorporating viscous drag effects. While it is tempting to invoke this effect for the differences shown in 
Figure 8, the independent total thrust calculations discussed in Section II are also off by similar amounts, 
indicating the discrepancies which appear in Figure 8 are more likely due to numerical inaccuracies rather 
than actual physical effects. This result points out an additional area in the model which needs improvement, 
and reiterates the need for independent checks to verify model predictions. 

Thrust Coefficients. The thrust for self- field MPD thrusters may be expressed 43 : 

T=bJ 2 (52) 

where b is the thrust coefficient, defined as: 



and 6 may vary between 1/4 and 3/4, depending upon where the current attaches to the cathode. Equation 
(52) provides a fairly accurate estimate of the electromagnetic thrust for coaxial MPD thrusters where b is 
well defined, but is inaccurate for flared thruster geometries and is unreliable when the thrust has a large 
electrothermal component. Figure 9 reproduces experimentally determined values of the thrust coefficient for 
the HSFAT as a function of J 2 /m, with model predictions overlayed on the experimental data. The code is 
able to accurately predict the thrust coefficient over the full range of J 2 fra values. Figure 10 plots the thrust 
coefficient versus current to bring out the dependence on mass flow; again, there is good agreement between 
measured and predicted values for all mass flow rates. As reported by Gilland 40 ’ 41 , the thrust coefficient 
approaches a value of LI — 1.2 x 10~ 7 N/A 2 as the electromagnetic thrust component begins to dominate 
the flow. Figure 11 displays the thrust coefficient versus J 2 /m using only the predicted electromagnetic 
component of thrust, plotted together with the thrust coefficients calculated using the total thrust. The 
electromagnetic thrust coefficients remain fairly constant at a value of approximately 1.1-1. 2xl0 -7 N/A 2 , 
indicating that the higher thrust coefficients observed at lower J 2 /ttl values are due primarily to large 
electrothermal thrust contributions rather than current redistribution along the flared anode or cathode 
surfaces. 

Visco us Lo sses. Figure 12 plots the ratio of calculated viscous drag to calculated total thrust as a function 
of J 2 /rh. The viscous drag ratio increases with increasing J 2 /m and decreasing mass flow rate, in qualitative 
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agreement with the scaling arguments of Wolff 44 . At a mass flow rate of 3 g/s, the viscous drag ratio increases 
from just over 1% of the total thrust at J 2 jrn = 2.5xl0 10 A 2 -s/kg to roughly 3.5% of the total thrust at 
1.5x 10 11 A 2 -s/kg. At the lower mass flow rate of 0.75 g/s, the viscous drag ratio rises from just over 2% 
of the total thrust at 2.5xl0 10 A 2 -s/kg to over 6% of the total thrust at 1.5xlO n A 2 -s/kg. These results 
indicate that some benefit in performance might be derived by operating at higher mass flow rates for a 
given J 2 jrn. 

Encl osed Current Contours. The enclosed current contours measured by Gilland 40 ’ 41 for the quasi-steady 
HSFAT are presented in Figure 13, for discharge currents of 7.9 kA and 17.8 kA and a mass flow rate of 
3 g/s. At 7.9 kA, roughly 80% of the current attaches along the straight anode section downstream of 
the anode flare, with 60% congregating along the last 2 cm of the anode. At the higher discharge current, 
nearly 40% of the current remains in the straight anode channel near the backplate, with the remaining 
60% spread along the anode downstream of the flared section. For comparison, the calculated enclosed 
current contours are presented as a matrix of mass flow rates and J 2 jrn values in Figure 14. There is 
little resemblance between the measured quasi-steady current distributions and the calculated steady-state 
current distributions. As discussed in Section Ilia, similar differences were reported by Sleziona et a/ 28,30 
when they compared calculated steady-state ZT-1 thruster current distributions with quasi-steady ZP-i 
thruster current distributions measured under different operating conditions but similar J 2 jrn . They argue 
that the current preferentially attaches near the cathode base in pulsed, quasi-steady thrusters, and is 
more evenly distributed along the hot cathode during steady-state thruster operation; consequently, the 
current distributions are not expected to be comparable. The use of different operating conditions tends to 
weaken this argument; however, both the Stuttgart code and the LeRC code predict significantly different 
steady-state current distributions than those actually observed in quasi-steady thruster operation, while 
accurately simulating anode current distributions in steady-state segmented anode thrusters. This raises the 
issue of whether quasi-steady current distribution measurements accurately represent steady-state current 
distributions, and the equally valid issue of whether numerical models which ignore electrode processes can 
accurately simulate current distributions within the thruster. 

The ability of the steady-state code to predict the quasi-steady HSFAT thrust and thrust coefficient val- 
ues under a variety of operating conditions indicates that the current distribution is less important than 
the total discharge current when calculating these global thruster parameters, and either quasi-steady or 
steady-state experimental data may be used to verify such predictions. The results of the ZT-1 thruster 
modeling with and without insulating segments, however, indicates that the current distribution does af- 
fect the development of the plasma flow properties (temperature, pressure, etc.), and underscores the need 
to determine whether quasi-steady thruster measurements are truly representative of steady-state thruster 
operation. An experimental comparison of anode current distributions and plasma plume properties for a 
single, segmented anode thruster operated under pulsed and steady-state conditions should be performed 
to expand upon the results of the Stuttgart ZT-l/ZP-1 current distribution experiments, and to determine 
the extent to which plasma flow conditions might vary between quasi-steady and steady-state operation. In 
addition, electrode processes must be coupled into the numerical model to evaluate the effect of sheaths on 
electrode current attachment. Little confidence can be placed in detailed comparisons between quasi-steady 
plasma flow measurements and numerical steady-state flow predictions until these issues are resolved. 

While thrust chamber predictions must await validation, useful trends can be established to provide insight 
into steady-state thruster operation. Figure 14 shows the predicted steady-state enclosed current contours 
for the HSFAT for a variety of mass flow rates and J 2 jrn values. The average magnetic Reynolds number 
for each case, defined by: 

Rm = iMjd-v cx L ( 54 ) 

where & is the average electrical conductivity, v cx is the exhaust velocity, and L is a length parameter, 
taken to be the difference between the anode and cathode radius, is displayed under each figure. Typically 
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40-50% of the current is contained within the upstream straight channel section of the thruster for all cases 
considered. For a given mass flow rate, a larger fraction of the discharge current is blown downstream of 
the anode exit plane as J 2 /m is increased. This result is commensurate with the increasing value of the 
magnetic Reynolds number, which is a measure of the ratio of magnetic convection to magnetic diffusion; as 
Rm increases, more of the magnetic field is convected with the plasma flow. This effect is observed in Figure 
14, where the enclosed current contours are equivalent to induced magnetic field streamlines (cf Section Ila). 
For a given J 2 /m, the magnetic Reynolds number decreases with decreasing mass flow rate, and the current 
is pulled back into the thruster. 

Figure 15 displays the fraction of electromagnetic thrust generated by current blown downstream of the 
anode exit plane for the HSFAT under the assumed steady-state operating conditions. The amount of 
electromagnetic thrust produced beyond the anode exit plane increases with increasing mass flow rate and 
increasing / 2 /m, in agreement with the predictions of Figure 14. At 3 g/s and 1.5 x 10 11 A 2 -s/kg, the amount 
of electromagnetic thrust created downstream of the anode is nearly 6% of the total electromagnetic thrust. 
Presumably this thrust is coupled through the magnetic field lines back into the thruster, but the details of 
how efficiently this process might occur remain to be explored. In addition, the downstream current may 
deposit thermal energy into the plume 45 rather than the chamber plasma, resulting in a loss of electrothermal 
thrust. 

Voltage and Efficiency. The radial electric field (Eqn. 7a) was numerically integrated from cathode to 
anode to provide an estimate of the potential drop across the plasma. The voltage calculations do not 
include fall voltages, but are still useful to predict voltage-current trends. Figure 16 shows the calculated 
plasma potential drop, and the experimentally determined total discharge voltage. The calculated plasma 
drop increases with increasing current, though at a slightly reduced rate compared to the total discharge 
voltage. In addition, the calculated potentials reproduce the dependence on mass flow rate seen in the total 
voltage measurements. 

As noted in the ZT-1 thruster comparisons, the calculated plasma potential depends strongly on the value 
used for the electrical conductivity. A reduced conductivity will raise the calculated value of the potential 
drop across the plasma, without strongly affecting the calculated thrust. To evaluate this effect under 
extreme variations, the electrical conductivity expression used in the code (Equation 3) was reduced by a 
factor of 10, and simulations were performed for J 2 /rh values of 5.1 x lO 10 and 1.5 x 10 11 A 2 -s/kg at a mass 
flow rate of 3 g/s. Thrust values increased a modest 20% in both cases, but the plasma potentials increased 
by factors of 4 to 5, bringing them close to the actual discharge voltages plotted in Figure 16. Changes in the 
plasma conductivity would thus appear to have a significant impact on total voltage, but minimal impact on 
the thrust. Choueiri et aP 6 have shown that current driven instabilities may be generated in MPD thrusters, 
with an attendant decrease in plasma conductivity. Under such conditions, the total discharge voltage would 
be increased without a commensurate increase in thrust. The thrust efficiency, defined 

VT = T 2 /(2mP) (55) 

where T is thrust, m is propellant mass flow rate, and P is power, would decrease under such operating 
conditions. The study of plasma instabilities and their effects on MPD thruster performance is clearly an 
area which warrants further investigation. 

Using the potential calculations displayed in Figure 16, the amount of power deposited into the plasma 
can be compared with the total thruster power, determined experimentally. Due to uncertainties in the 
electrical conductivity a direct comparison is of limited use, but trends may be discerned which aid in 
thruster performance evaluation. Figure 17 displays the amount of power deposited into the plasma and the 
total power deposited into the thruster as a function of J 2 /m. The predicted trends in power with increasing 
J 2 /m are similar for both the total thruster power and the power deposited solely into the plasma for all mass 
flow rates considered. Using the expression for efficiency defined above, both the total thruster efficiency (*7 t ) 
and the efficiency with which power deposited in the plasma is converted to thrust (frozen flow efficiency, 
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r\ff), can be separately evaluated. Figure 18 displays the experimentally evaluated thruster efficiency and 
the calculated frozen flow efficiency. The total thruster efficiency stays fairly flat in the 5-15% range, with 
the higher mass flow rates providing slightly higher thruster efficiencies. The frozen flow efficiencies are 
highest at the lower J 2 /rn values, ranging from 40% at m=1.5 g/s to 60% at m=3 g/s. The frozen flow 
efficiencies decline with increasing J 2 /m, leveling out at a value around 37% at the higher J 2 jm values. 
Assuming a smaller value for the electrical conductivity would result in even lower frozen flow efficiencies, 
because the amount of power deposited in the plasma would increase without significantly increasing the 
thrust. These results, while subject to the numerous approximations made in the model, predict that low 
frozen flow efficiencies will limit self-field MPD thruster performance. The impact of thruster geometry and 
propellant species on the frozen flow efficiency is under investigation. 


IV. Concluding Remarks, 


A two dimensional magnetohydrodyanamics code has been developed to predict steady-state MPD thruster 
performance. The governing equations and corresponding methods of solution were outlined and discussed. 
Experimental comparisons were used to evaluate model calculations for two thruster geometries over a range 
of discharge currents and propellant mass flow rates. The self-field model accurately predicted thruster 
performance up to J 2 /m values slightly below the onset of measured voltage instabilities. Methods to 
improve code performance have been outlined and discussed. Nonuniform grids should be incorporated to 
better model the sharp spatial gradients in temperature and pressure near the electrode surfaces. Improved 
estimates for the plasma transport coefficients and accurate electrode sheath models must be incorporated to 
accurately predict plasma potentials and total discharge voltages. Issues concerning the use of quasi-steady 
thruster data to verify steady-state model predictions must be addressed. 
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Figure 1 . Self- field MPD thruster diagram. 



Figure 2. MPD thruster program outline. 
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Figure 3(a). ZT-1 thruster schematic 



Figure 3(b). U. Stuttgart enclosed current predictions. Contour levels every 10%. 
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Figure 3(c). NASA LeRC enclosed current predictions. Contour levels every 10%. 


Figure 3. Predicted enclosed current contour comparisons for the ZT-1 thruster operat- 
ing at 6000 A, 6 g/s, argon propellant. 
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Figure 4. Predicted (a) temperature, (b) pressure, and (c) density contours for ZT-1 
thruster at 6000 A, 6 g/s, argon propellant. Pressure and density are normalized to inlet 
values: P 0 = l.l x 10 3 N/m 2 , p o = 1.8xl0" 3 kg/m 3 . Plots (b) and (c) are logio of normalized 
values. 
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Figure 4, cont. Predicted (d) electric potential, (e) Mach number, and (f) velocity flow 
distribution for ZT-1 thruster operating at 6000 A, 6 g/s. 
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Figure 5. Predicted (a) temperature, (b) pressure, and (c) density contours for ZT-1 
thruster without insulated anode segments, operating at 6000 A, 6 g/s, argon propellant. 
Pressure and density are normalized to inlet values: P o =l.lxl0 3 N/m 2 , p 0 ~l^Sx 10 -3 kg/m 3 . 
Plots (b) and (c) are logio of normalized values. 
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Figure 5, cont. Predicted (d) electric potential, (e) Mach number, and (f) velocity flow 
distribution for ZT-l thruster without insulated anode segments, operating at 6000 A, 6 g/s. 
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Figure 5(g) Predicted enclosed current contours for the ZT-1 thruster without insulating 
anode segments, operating at 6000 A and 6 g/s, argon propellant. 



Figure 6. Half-scale flared anode thruster (HSFAT) diagram. 
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Figure 8. Calculated electromagnetic and total thrust for HSFAT, argon propellant. 
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Figure 10. HSFAT thrust coefficients plotted against current to show mass flow rate 
dependence. 
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Figure 12. Ratio of viscous drag to total thrust, HSFAT thruster, argon propellant. 
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Figure 13 . Measured enclosed current contours, HSFAT, argon propellant. 
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Figure 15. Predicted precentage of electromagnetic thrust generated downstream of the 
HSFAT anode. 



Figure 16- Measured total voltage and predicted plasma voltage for HSFAT thruster, 
argon propellant. 
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Figure 18 . Measured total thruster efficiency and predicted frozen flow efficiency (plasma) 
for HSFAT thruster, argon propellant. 
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